#' bngeneplot
#'
#' Plot gene relationship within the specified pathway
#'
#' @param results the enrichment analysis result from clusterProfiler or ReactomePA
#' @param exp gene expression matrix
#' @param expRow the type of the identifier of rows of expression matrix
#' @param expSample candidate samples to be included in the inference, default to all
#' @param algo structure learning method used in boot.strength(), default to "hc"
#' @param algorithm.args parameters to pass to bnlearn structure learnng function
#' @param R the number of bootstrap
#' @param pathNum the pathway number (the number of row of the original result, ordered by p-value)
#' @param convertSymbol whether the label of resulting network is converted to symbol, default to TRUE
#' @param interactive whether to use bnviewer (default to FALSE)
#' @param cexCategory scaling factor of size of nodes
#' @param delZeroDegree delete zero degree nodes
#' @param disc discretize the expressoin data
#' @param tr Specify data.frame if one needs to discretize as the same parameters as the other dataset
#' @param remainCont Specify characters when perform discretization, if some columns are to be remain continuous
#' @param cl cluster object from parallel::makeCluster()
#' @param showDir show the confidence of direction of edges
#' @param showDepHist whether to show depmap histogram
#' @param chooseDir if undirected edges are present, choose direction of edges (default: FALSE)
#' @param scoreType score type to use on choosing direction
#' @param labelSize the size of label of the nodes
#' @param layout ggraph layout, default to "nicely"
#' @param clusterAlpha if specified multiple pathways, the parameter is passed to geom_mark_hull()
#' @param strType "normal" or "ms" for multiscale implementation
#' @param sp query to graphite::pathways(), default to "hsapiens"
#' @param compareRef whether compare to the reference network
#' @param compareRefType "intersection" or "difference"
#' @param pathDb query to graphite::pathways(), default to "reactome"
#' @param dep the tibble storing dependency score from library depmap
#' @param sizeDep whether to reflect DepMap score to the node size
#' @param cellLineName the cell line name to be included
#' @param strengthPlot append the barplot depicting edges with high strength
#' @param nStrength specify how many edges are included in the strength plot
#' @param showLineage show the dependency score across the lineage
#' @param depMeta depmap::depmap_metadata(), needed for showLineage
#' @param strThresh the threshold for strength
#' @param hub visualize the genes with top-n hub scores
#' @param returnNet whether to return the network
#' @param otherVar other variables to be included in the inference
#' @param otherVarName the names of other variables
#' @param onlyDf return only data.frame used for inference
#' @param orgDb perform clusterProfiler::setReadable based on this organism database
#' @param shadowText whether to use shadow text for the better readability (default: TRUE)
#' @param bgColor color for text background when shadowText is TRUE
#' @param textColor color for text when shadowText is TRUE
#' @param seed A random seed to make the analysis reproducible, default is 1.
#' @return ggplot2 object
#'
#' @examples
#' data("exampleEaRes");data("exampleGeneExp")
#' res <- bngeneplot(results = exampleEaRes, exp = exampleGeneExp, pathNum = 1, R = 10, convertSymbol = TRUE, expRow = "ENSEMBL")
#'
#' @importFrom dplyr group_by summarize arrange n
#' @importFrom graphite pathways convertIdentifiers pathwayGraph
#' @importFrom clusterProfiler setReadable
#' @export
#'

bngeneplot <- function (results, exp, expSample=NULL, algo="hc", R=20, returnNet=FALSE, algorithm.args=NULL,
                        pathNum=NULL, convertSymbol=TRUE, expRow="ENSEMBL", interactive=FALSE,
                        cexCategory=1, cl=NULL, showDir=FALSE, chooseDir=FALSE, scoreType="bic-g",
                        labelSize=4, layout="nicely", clusterAlpha=0.2, strType="normal", delZeroDegree=TRUE,
                        otherVar=NULL, otherVarName=NULL, onlyDf=FALSE, disc=FALSE, tr=NULL, remainCont=NULL,
                        sp="hsapiens", compareRef=FALSE, compareRefType="intersection", pathDb="reactome",
                        dep=NULL, depMeta=NULL, sizeDep=FALSE, showDepHist=TRUE, cellLineName="5637_URINARY_TRACT",
                        showLineage=FALSE, orgDb=org.Hs.eg.db, shadowText=TRUE, bgColor="white", textColor="black",
                        strengthPlot=FALSE, nStrength=10, strThresh=NULL, hub=NULL, seed = 1) {
    
    if (is.null(expSample)) {expSample=colnames(exp)}
    if (compareRef & length(pathNum) > 1){stop("compareRef can be used with one pathNum or pathName.")}
    if (interactive & compareRef){stop("compareRef must be set to FALSE when use bnviewer.")}
    if (!is.numeric(pathNum)){stop("Please specify number(s) for pathNum.")}
    if (sizeDep & !convertSymbol){stop("sizeDep must be used with convertSymbol set to TRUE.")}
    if (sizeDep & length(pathNum) > 1){stop("sizeDep can be used with one pathNum.")}
    if (compareRef & is.null(pathDb)){stop("please specify which database to use as reference.")}
    if (showLineage & strengthPlot){stop("please specify one of showLineage or strengthPlot.")}
    if (results@ontology=="Reactome") {
        if (grepl("Homo sapiens\r: ", results@result$Description[1], fixed = TRUE))
            results@result$Description <- vapply(strsplit(results@result$Description, "\r: "),
                "[", 2, FUN.VALUE="character")
    }
    # if (results@keytype == "kegg"){
    #     resultsGeneType <- "ENTREZID"
    # } else {
    #     resultsGeneType <- results@keytype
    # }
    if (!is.null(orgDb)){
        results <- setReadable(results, OrgDb=orgDb)
    }
    tmpCol <- colnames(results@result)
    tmpCol[tmpCol=="core_enrichment"] <- "geneID"
    tmpCol[tmpCol=="qvalues"] <- "qvalue"
    tmpCol[tmpCol=="setSize"] <- "Count"
    colnames(results@result) <- tmpCol

    if (showLineage) {
        if (is.null(dep)){dep = depmap::depmap_crispr()}
        if (is.null(depMeta)){depMeta = depmap::depmap_metadata()}
    }

    if (sizeDep) {
        if (is.null(cellLineName)){stop("Please specify cell line name.")}
        if (is.null(dep)){dep = depmap::depmap_crispr()}
        filteredDep <- dep %>% filter(cell_line==cellLineName)
        depHist <- ggplot(filteredDep, aes(x=dependency)) +
            geom_histogram(aes(fill=..count..), col="black") +
            scale_fill_gradient("Count", low = "blue", high = "red") +
            theme_minimal(base_family = "Arial Narrow") +
            ggtitle(cellLineName)+
            theme(plot.title = element_text(hjust=0.5, face="bold"),
                  axis.text = element_text(size=10),
                  axis.title = element_text(size=12))
    }
    if (sizeDep){
        ## when reflecting DepMap scores to node size
        scaleSizeLow <- 1
        scaleSizeHigh <- 10
    } else {
        scaleSizeLow <- 3
        scaleSizeHigh <- 8
    }

    res <- results@result

    genesInPathway <- unique(unlist(strsplit(res[pathNum, ]$geneID, "/")))
    if (!is.null(orgDb)) {
        genesInPathway <- clusterProfiler::bitr(genesInPathway,
                                                fromType="SYMBOL",
                                                toType=expRow,
                                                OrgDb= orgDb )[expRow][,1]
    }
    pcs <- exp[ intersect(rownames(exp), genesInPathway), expSample ]

    if (expRow!="SYMBOL"){
        if (convertSymbol && !is.null(orgDb)) {
          # rn <- clusterProfiler::bitr(rownames(pcs),
          #                             fromType=expRow,
          #                             toType="SYMBOL",
          #                             OrgDb=org.Hs.eg.db)["SYMBOL"][,1]
          ## Change expression matrix rownames to symbol
          ## If one "expRow" hit to multiple symbols, delete the ID from the subsequent analysis, showing warning.
          matchTable <- clusterProfiler::bitr(rownames(pcs), fromType=expRow, toType="SYMBOL", OrgDb=orgDb)
          if (sum(duplicated(matchTable[,1])) >= 1) {
            message("Removing IDs that matches the multiple symbols")
            matchTable <- matchTable[!matchTable[,1] %in% matchTable[,1][duplicated(matchTable[,1])],]
          }
          rnSym <- matchTable["SYMBOL"][,1]
          rnExp <- matchTable[expRow][,1]
          pcs <- pcs[rnExp, ]
          rownames(pcs) <- rnSym
        }
    }

    pcs <- data.frame(t(pcs))

    if (sizeDep){
        pcs <- pcs[,intersect(filteredDep$gene_name, colnames(pcs)),]
    }


    if (length(pathNum) >= 2) {
        clus <- c()

        for (num in pathNum){
            tmp <- data.frame(strsplit(res[num, ]$geneID, "/")[[1]])
            tmp$Pathway <- res[num, ]$Description
            clus <- rbind(clus, tmp)
        }

        colnames(clus) <- c("geneID", "Pathway")

        if (convertSymbol) {
            # m <- clusterProfiler::bitr(clus$geneID,
            #                            fromType=resultsGeneType,
            #                            toType="SYMBOL",
            #                            OrgDb=org.Hs.eg.db)
            # colnames(m) <- c("geneID","SYMBOL")
            mcl <- clus#merge(m, clus)
            colnames(mcl) <- c("SYMBOL","Pathway")
            cnt <- mcl %>% group_by(SYMBOL) %>%
                arrange(Pathway, .by_group = TRUE) %>%
                summarize(n=n(), Pathway=paste0(Pathway, collapse = " + "))
            ovl <- cnt[cnt$n > 1,]
            mclSub <- subset(mcl, mcl$SYMBOL %in% cnt[cnt$n==1,]$SYMBOL)
            cls <- rbind(mclSub[,c("SYMBOL","Pathway")], ovl[,c("SYMBOL","Pathway")])
            rownames(cls) <- cls$SYMBOL

        } else {
            if (!is.null(orgDb)){
                m <- clusterProfiler::bitr(clus$geneID,
                                           fromType="SYMBOL",
                                           toType=expRow,
                                           OrgDb=orgDb)
                colnames(m) <- c("geneID", expRow)
                mcl <- merge(m, clus)
                cnt <- mcl %>% group_by_at(expRow) %>%
                    arrange(Pathway, .by_group = TRUE) %>%
                    summarize(n=n(), Pathway=paste0(Pathway, collapse = " + "))
                ovl <- cnt[cnt$n > 1,]
                mclSub <- subset(mcl, mcl[expRow][,1] %in% as.character(as.matrix(cnt[cnt$n==1,][,expRow])))
                cls <- data.frame(rbind(mclSub[,c(expRow,"Pathway")], ovl[,c(expRow,"Pathway")]))
                rownames(cls) <- cls[, 1]
            }
        }
    }
    
    geneNames <- colnames(pcs)

    ## Insert other vars
    if (!is.null(otherVar)) {
        pcs <- cbind(pcs, otherVar)
        if (!is.null(otherVarName)){
          colnames(pcs) <- c(geneNames, otherVarName)
        }
    }

    if (disc){
        pcsRaw <- pcs
        pcs <- discDF(pcs, tr=tr, remainCont=remainCont)
    } else {
        pcsRaw <- pcs ## Hold pcsRaw anyway
    }

    if (onlyDf){
        return(pcs)
    }

    # print(dim(pcs))
    if (dim(pcs)[2]<=1){message("the number of gene is zero or one");return("error")}

    ## Bootstrap-based inference
    if (strType == "normal"){
      strength <- withr::with_seed(seed = seed, boot.strength(pcs, algorithm=algo, algorithm.args=algorithm.args, R=R, cluster=cl))
    } else if (strType == "ms"){
      strength <- withr::with_seed(seed = seed, inferMS(pcs, algo=algo, algorithm.args=algorithm.args, R=R, cl=cl))
    }

    ## Barplot of edge strength
    if (strengthPlot){
        strengthTop <- strength[order(strength$strength+strength$direction, decreasing = TRUE),][seq_len(nStrength),]
        strengthTop$label <- paste(strengthTop$from, "->", strengthTop$to)
        stp <- strengthTop %>% tidyr::pivot_longer(cols=c(strength, direction)) %>%
            ggplot(aes(x=label, y=value, fill=name))+geom_bar(position="dodge",stat="identity")+
            coord_flip(y=c(min(strengthTop$strength, strengthTop$direction)-0.05,1.0))+xlab("edges")+
            theme_bw()+scale_fill_manual(values = c("tomato","dodgerblue")) +
            scale_x_discrete(labels = function(x) stringr::str_wrap(x, width = 25))
    }

    ## Average by specified threshold
    if (!is.null(strThresh)){
        av <- averaged.network(strength, threshold=strThresh)
    } else {
        av <- averaged.network(strength)
    }

    # if (chooseDir){
    #     av <- chooseEdgeDir(av, pcs, scoreType)
    # }

    av <- cextend(av, strict=FALSE)

    if (interactive) {
        strength.viewer(bayesianNetwork = av, bayesianNetwork.boot.strength = strength)
    } else {

        g <- bnlearn::as.igraph(av)
        e <- as_edgelist(g, names = TRUE)
        if (dim(e)[1]==0){message("no edge present in graph");return("error")}
        eName <-paste0(e[,1], "_", e[,2])
        colnames(e) <- c("from","to")
        eDf <- merge(e, strength)
        rownames(eDf) <- paste0(eDf$from, "_", eDf$to)
        eDf <- eDf[eName, ]
        g <- set.edge.attribute(g, "color", index=E(g), eDf$strength)
        if (showDir){
            g <- set.edge.attribute(g, "label", index=E(g), round(eDf$direction,2))
        } else {
            g <- set.edge.attribute(g, "label", index=E(g), NA)
        }
        
        ## Hub genes
        hScore <- hub.score(g, scale = TRUE, weights = E(g)$color)$vector

        if (!is.null(hub)){
            defHub <- hScore[order(hScore, decreasing=TRUE)][seq_len(hub)]
            nodeShape <- names(V(g)) %in% names(defHub)
            nodeShape <- ifelse(nodeShape, 19, 21)
            V(g)$shape <- nodeShape
        } else {
            V(g)$shape <- rep(19, length(V(g)))
        }

        E(g)$width <- E(g)$color
        edgeWName <- "strength"

        if (sizeDep){
            sizeLab <- "-dependency"
            filteredDep <- filteredDep %>% filter(gene_name %in% names(V(g))) %>%
                arrange(match(gene_name, names(V(g))))
            depHistSub <- ggplot(filteredDep, aes(x=dependency)) +
                geom_histogram(binwidth=0.5, aes(fill=..count..), col="black") +
                scale_fill_gradient("Count", low = "blue", high = "red")+
                theme_minimal(base_family = "Arial Narrow") +
                theme(axis.text=element_text(size=10), axis.title=element_text(size=12))
            # Subset to those dependency scores are available
            depSubG <- V(g)[names(V(g)) %in% c(filteredDep$gene_name) | names(V(g)) %in% tail(colnames(pcs), n=dim(otherVar)[2])]
            g <- igraph::subgraph(g, depSubG)
            tmpSize <- vapply(names(V(g)), function(x) ifelse(x %in% filteredDep$gene_name, -1 * as.numeric(subset(filteredDep, gene_name==x)$dependency), NA), FUN.VALUE=1) #-1 * filteredDep$dependency
            # Size of metadata is set to mean of score
            tmpSize[is.na(tmpSize)] <- mean(tmpSize, na.rm=TRUE)
            V(g)$size <- tmpSize
            #meanExp <- apply(pcs[, names(V(g))], 2, mean)
            meanExpCol <- vapply(names(V(g)), function(x) ifelse(x %in% geneNames, mean(pcsRaw[, x]), NA), FUN.VALUE=1)
        } else {
            sizeLab <- "expression"
            #meanExp <- apply(pcs[, names(V(g))], 2, mean)
            meanExpCol <- vapply(names(V(g)), function(x) ifelse(x %in% geneNames, mean(pcsRaw[, x]), NA), FUN.VALUE=1)
            meanExpSize <- meanExpCol
            meanExpSize[is.na(meanExpSize)] <- 1
            V(g)$size <- meanExpSize
        }

        ## Node color
        V(g)$color <- meanExpCol

        ## Cluster for multiple pathways
        if (length(pathNum) > 1) {
            V(g)$Pathway <- vapply(names(V(g)), function(x) ifelse(x %in% geneNames, cls[x, ]$Pathway, "other variables"), FUN.VALUE="character")
        }

        ## Plot
        if (length(pathNum) == 1) {
            if (delZeroDegree){
                delG <- delete.vertices(g, igraph::degree(g)==0)
            } else {
                delG <- g
            }
            p <- ggraph(delG, layout=layout) + geom_edge_diagonal(edge_alpha=1,
                                                               position="identity",
                                                               aes_(edge_colour=~I(color), width=~I(width), label=~I(label)),
                                                               label_size=3*(labelSize/4),
                                                               label_colour=NA,
                                                               angle_calc = "along",
                                                               label_dodge=unit(3,'mm'),
                                                               arrow=arrow(length=unit(4, 'mm')),
                                                               end_cap=circle(5, 'mm'))+
                geom_node_point(aes_(color=~color, size=~size, shape=~shape), show.legend=TRUE)+
                scale_color_continuous(low="blue", high="red", name="expression") +
                scale_size(range=c(scaleSizeLow, scaleSizeHigh) * cexCategory, name=sizeLab)+
                scale_edge_width(range=c(1, 3), guide="none")+
                scale_edge_color_continuous(low="dodgerblue", high="tomato", name="strength")+
                guides(edge_color = guide_edge_colorbar(title.vjust = 3))+
                # geom_node_text(aes_(label=~name), check_overlap=TRUE, repel=TRUE, size = labelSize) +
                scale_shape_identity()+
                theme_graph() + ggtitle(res[pathNum, "Description"])

            if (shadowText){
                p <- p + geom_node_text(aes_(label=~stringr::str_wrap(name, width = 25)),
                    check_overlap=TRUE, repel=TRUE, size = labelSize,
                    color = textColor,
                    bg.color = bgColor, segment.color="black",
                    bg.r = .15)
            } else {
                p <- p + geom_node_text(aes_(label=~stringr::str_wrap(name, width = 25)),
                    check_overlap=TRUE, repel=TRUE, size = labelSize)
            }

            if (sizeDep & !compareRef){
                if (showDepHist){
                    layoutSizedep <- "
                          ACC
                          BCC
                          "
                    p <- depHist + depHistSub + p + plot_layout(design=layoutSizedep)
                }
            }

            if (compareRef){
                pathName = res[pathNum, "Description"]
                graphiteP = pathways(sp, pathDb)[[pathName]]
                if (!convertSymbol) {
                    graphiteP <- convertIdentifiers(graphiteP, expRow)
                } else {
                    graphiteP <- convertIdentifiers(graphiteP, "symbol")
                }

                refG <- pathwayGraph(graphiteP)
                refG <- igraph.from.graphNEL(refG, name=TRUE)
                refG <- set.vertex.attribute(refG, "name", value=vapply(strsplit(names(V(refG)), ":"), "[", 2, FUN.VALUE=character(1)))

                refV <- names(V(refG))
                intNodeLen <- length(intersect(refV, names(V(g))))

                refVSubset <- V(refG)[intersect( refV, names(V(g)))]
                avVSubset <- V(g)[intersect( refV, names(V(g)))]

                refSubG <- igraph::subgraph(refG, refVSubset)
                avSubG <- igraph::subgraph(g, avVSubset)

                refSubGUnd <- as.undirected(refSubG)
                avSubGUnd <- as.undirected(avSubG)

                difG <- difference(g, refSubG)
                difG <- delete.vertices(difG, igraph::degree(difG)==0)

                refELen <- length(E(refSubG))
                avELen <- length(E(avSubG))

                intG <- intersection(refSubG, avSubG, keep.all.vertices = FALSE)
                intG <- delete.vertices(intG, igraph::degree(intG)==0)


                if (compareRefType == "intersection"){
                    refPlot <- intG
                    t <- "Overlapping"
                } else if (compareRefType == "difference"){
                    refPlot <- difG
                    t <- "Different"
                }

                ovlELen <- length(E(refPlot))
                intP <- ggraph(refPlot, layout=layout) + geom_edge_diagonal(edge_alpha=1,
                                                                            position="identity",
                                                                            aes_(edge_colour=~I(color), width=~I(width), label=~I(label)),
                                                                            label_size=3*(labelSize/4),
                                                                            label_colour=NA,
                                                                            angle_calc = "along",
                                                                            label_dodge=unit(3,'mm'),
                                                                            arrow=arrow(length=unit(4, 'mm')),
                                                                            end_cap=circle(5, 'mm'))+
                    geom_node_point(aes_(color=~color, size=~size, shape=~shape), show.legend=TRUE)+
                    scale_color_continuous(low="blue", high="red", name = "expression") +
                    scale_size(range=c(scaleSizeLow, scaleSizeHigh) * cexCategory, name=sizeLab)+
                    scale_edge_width(range=c(1, 3), guide="none")+
                    scale_edge_color_continuous(low="dodgerblue", high="tomato", name="strength")+
                    guides(edge_color = guide_edge_colorbar(title.vjust = 3))+
                    # geom_node_text(aes_(label=~name), check_overlap=TRUE, repel=TRUE, size = labelSize) +
                    theme_graph() +
                    scale_shape_identity()+
                    ggtitle(paste(paste("Overlapping nodes =", intNodeLen),
                                  paste("Edge number in reference network =", refELen),
                                  paste("Edge number in inferred network =", avELen),
                                  paste(paste0(t, " edges ="), ovlELen),
                                  sep="\n"))
                if (shadowText){
                    intP <- intP + geom_node_text(aes_(label=~stringr::str_wrap(name, width = 25)),
                        check_overlap=TRUE, repel=TRUE, size = labelSize,
                        color = "white",
                        bg.color = "black", segment.color="black",
                        bg.r = .15)
                } else {
                    intP <- intP + geom_node_text(aes_(label=~stringr::str_wrap(name, width = 25)),
                        check_overlap=TRUE, repel=TRUE, size = labelSize)
                }

                p2 <- p + theme(legend.position="none")
                if (sizeDep & showDepHist){
                    layoutDep <- "
                          ACCDD
                          BCCDD
                          "
                    p <- depHist + depHistSub + p2 + intP + plot_layout(design=layoutDep)

                } else {
                    p <- p2 + intP
                }
            }

            if (showLineage){
                lineageP = depMeta %>%
                    dplyr::select(depmap_id, lineage) %>%
                    dplyr::full_join(dep, by = "depmap_id") %>%
                    dplyr::filter(gene_name %in% names(V(g))) %>%
                    ggplot(aes(x=lineage, y=dependency, fill=lineage)) +
                    geom_boxplot(outlier.alpha = 0.1) +
                    theme_bw()+
                    theme(axis.text.x = element_text(angle = 45, hjust=1),
                          axis.text = element_text(size=12),
                          axis.title = element_text(size=14),
                          legend.position = "none")
                p <- p / lineageP

            }
        } else if (length(pathNum) > 1) {
            if (delZeroDegree){
                delG <- delete.vertices(g, igraph::degree(g)==0)
            } else {
                delG <- g
            }
            xy <- graphlayouts::layout_as_backbone(igraph::as.undirected(delG))$xy
            p <- ggraph(delG, layout="manual", x=xy[,1], y=xy[,2]) + geom_edge_diagonal(edge_alpha=1,
                                                                                     position="identity",
                                                                                     aes_(edge_colour=~I(color), width=~I(width), label=~I(label)),
                                                                                     label_size=3*(labelSize/4),
                                                                                     label_colour=NA,
                                                                                     angle_calc = "along",
                                                                                     label_dodge=unit(3,'mm'),
                                                                                     arrow=arrow(length=unit(4, 'mm')),
                                                                                     end_cap=circle(5, 'mm'))+
                geom_node_point(aes_(color=~color, size=~size, shape=~shape), show.legend=TRUE)+
                ggforce::geom_mark_hull(
                    aes(xy[,1], xy[,2], group = Pathway, fill = Pathway, label=Pathway),
                    concavity = 8, expand = unit(2, "mm"), alpha = clusterAlpha, label.fill="transparent",
                    show.legend = FALSE, label.fontsize=12) +
                scale_fill_discrete(guide="none")+
                scale_color_continuous(low="blue", high="red", name = "expression") +
                scale_size(range=c(3, 8) * cexCategory, name=sizeLab)+
                scale_edge_width(range=c(1, 3), guide="none")+
                scale_edge_color_continuous(low="dodgerblue", high="tomato", name="strength")+
                guides(edge_color = guide_edge_colorbar(title.vjust = 3))+
                scale_shape_identity()+
                geom_node_text(aes_(label=~name), check_overlap=TRUE, repel=TRUE, size = labelSize) +
                theme_graph()
        }

        if (strengthPlot){
            p <- p / stp + plot_layout(nrow=2, ncol=1, heights=c(0.8, 0.2))
        }

        if (returnNet){
            returnList <- list()
            returnList[["plot"]] <- p
            returnList[["str"]] <- strength
            returnList[["av"]] <- av
            returnList[["df"]] <- pcs
            return(returnList)
        } else {
            return(p)
        }
    }
}
